{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "831ce106-fe2f-45d4-a12b-50709043948e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9221f194-8590-4689-a4d7-279a27ab98c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scikit-learn: 1.2.2\n",
      "mapie       : 0.6.4\n",
      "numpy       : 1.23.5\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%watermark -p scikit-learn,mapie,numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d641978b-0ee5-4bf5-849c-5b575e1de6a9",
   "metadata": {},
   "source": [
    "## Define a Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9750c18e-807d-4615-8055-d4c52fc2f94a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "\n",
    "iris = load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "\n",
    "X_, X_test, y_, y_test = train_test_split(\n",
    "    X, y, test_size=0.015, random_state=123, stratify=y)\n",
    "    \n",
    "# make training and calibration dataset\n",
    "X_train, X_calib, y_train, y_calib = train_test_split(\n",
    "    X_, y_, test_size=0.5, random_state=123, stratify=y_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7426fced-f48b-4ca0-a535-a0d2d502391c",
   "metadata": {},
   "source": [
    "## Train a Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "585a361b-5c4b-482b-b1dc-715a51458004",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "model = LogisticRegression(random_state=123, solver=\"newton-cg\")\n",
    "\n",
    "model.fit(X_train, y_train);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4be80e-e960-4aaf-ad7d-b552196cec87",
   "metadata": {},
   "source": [
    "## Get Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "135e10e1-b701-4371-b8b8-c6e7e15fe10a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(74, 3)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = model.predict_proba(X_calib)\n",
    "\n",
    "# predictions have the shape [n_samples, n_classes]\n",
    "pred.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9ddad3d5-fd48-4cba-8679-8047d85e2d08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.21511883e-06, 2.02140907e-02, 9.79782694e-01])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# each entry has 3 probability scores here, 1 for each class\n",
    "pred[0] # first prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dc0bd0b-ccc7-4411-b43b-43bdd94e9c85",
   "metadata": {},
   "source": [
    "# Conformal Prediction From Scratch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc2108d7-c914-42d8-8f11-cfc73e0fb32b",
   "metadata": {},
   "source": [
    "#### Compute quantile threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "477c9253-ab5b-4cbc-92c1-a15c3b9041ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "confidence_level = 0.95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fa3d3f4c-6862-4a45-aae4-f8d8b4e334d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(74,)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Select the probabilities corresponding to the true class labels\n",
    "\n",
    "proba_true = pred[np.arange(y_calib.shape[0]), y_calib]\n",
    "print(proba_true.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2a2a56b9-a40c-41cd-8639-f8e5783ab6f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute nonconformity score\n",
    "\n",
    "nc_scores = 1 - proba_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cc46575a-8290-45f9-90b2-b8bbfd82450c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAteElEQVR4nO3de1TVdb7/8deWyxYRKDW5BKEmWF7LnLyVoimljZN6pqms1LROHq000o6OndEuR0yXjJWXnKZQZzJtTJ1m5Y0y8ZYpKtWooxxFwYRIM0BNTPn8/mi5f+1AhC2w94d5Ptb6rtX38/3sz/e9P4D71ef73Xs7jDFGAAAAlqrn7QIAAACuBmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBq/t4uoKaVlpbq+PHjCgkJkcPh8HY5AACgEowxKi4uVlRUlOrVq3jtpc6HmePHjysmJsbbZQAAAA/k5uYqOjq6wj51PsyEhIRI+mkyQkNDvVwNAACojKKiIsXExLhexytS58PMpUtLoaGhhBkAACxTmVtEuAEYAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDV/bxdgu2YTP6qRcY9Mv7dGxgUAoK5hZQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFbzapiZP3++2rdvr9DQUIWGhqpr165as2aN67gxRlOnTlVUVJSCgoKUkJCgvXv3erFiAADga7waZqKjozV9+nRlZGQoIyNDvXv31n333ecKLDNmzFBKSormzJmjnTt3KiIiQn379lVxcbE3ywYAAD7Eq2FmwIAB6t+/v+Lj4xUfH6///d//VcOGDbV9+3YZYzR79mxNnjxZgwcPVtu2bbVo0SKdPXtWS5Ys8WbZAADAh/jMPTMXL17U0qVLdebMGXXt2lXZ2dnKz89XYmKiq4/T6VTPnj21bdu2y45TUlKioqIitw0AANRdXg8zX331lRo2bCin06lRo0Zp5cqVat26tfLz8yVJ4eHhbv3Dw8Ndx8qTnJyssLAw1xYTE1Oj9QMAAO/yephp1aqVMjMztX37dv3Xf/2Xhg0bpn379rmOOxwOt/7GmDJtPzdp0iQVFha6ttzc3BqrHQAAeJ+/twsIDAxUy5YtJUmdOnXSzp079dprr+m///u/JUn5+fmKjIx09S8oKCizWvNzTqdTTqezZosGAAA+w+srM79kjFFJSYmaN2+uiIgIpaWluY6dP39e6enp6tatmxcrBAAAvsSrKzO///3v1a9fP8XExKi4uFhLly7Vxo0btXbtWjkcDo0bN07Tpk1TXFyc4uLiNG3aNDVo0EBDhgzxZtkAAMCHeDXMfPPNN3r00UeVl5ensLAwtW/fXmvXrlXfvn0lSc8//7x++OEHjR49WqdOnVLnzp21fv16hYSEeLNsAADgQxzGGOPtImpSUVGRwsLCVFhYqNDQ0Gofv9nEj6p9TEk6Mv3eGhkXAAAbVOX12+fumQEAAKgKwgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwmlfDTHJysn71q18pJCRETZs21cCBA3XgwAG3PsOHD5fD4XDbunTp4qWKAQCAr/FqmElPT9eYMWO0fft2paWl6cKFC0pMTNSZM2fc+t1zzz3Ky8tzbatXr/ZSxQAAwNf4e/Pka9euddtPTU1V06ZNtWvXLvXo0cPV7nQ6FRERUdvlAQAAC/jUPTOFhYWSpEaNGrm1b9y4UU2bNlV8fLyeeOIJFRQUXHaMkpISFRUVuW0AAKDu8pkwY4xRUlKS7rjjDrVt29bV3q9fP7377rvasGGDZs2apZ07d6p3794qKSkpd5zk5GSFhYW5tpiYmNp6CgAAwAscxhjj7SIkacyYMfroo4+0ZcsWRUdHX7ZfXl6eYmNjtXTpUg0ePLjM8ZKSEregU1RUpJiYGBUWFio0NLTa62428aNqH1OSjky/t0bGBQDABkVFRQoLC6vU67dX75m55Omnn9aHH36oTZs2VRhkJCkyMlKxsbHKysoq97jT6ZTT6ayJMgEAgA/yapgxxujpp5/WypUrtXHjRjVv3vyKjzl58qRyc3MVGRlZCxUCAABf59V7ZsaMGaO//vWvWrJkiUJCQpSfn6/8/Hz98MMPkqTTp09r/Pjx+uyzz3TkyBFt3LhRAwYMUJMmTTRo0CBvlg4AAHyEV1dm5s+fL0lKSEhwa09NTdXw4cPl5+enr776SosXL9b333+vyMhI9erVS8uWLVNISIgXKgYAAL7G65eZKhIUFKR169bVUjUAAMBGPvPWbAAAAE8QZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqHoWZ7Ozsajl5cnKyfvWrXykkJERNmzbVwIEDdeDAAbc+xhhNnTpVUVFRCgoKUkJCgvbu3Vst5wcAAPbzKMy0bNlSvXr10l//+ledO3fO45Onp6drzJgx2r59u9LS0nThwgUlJibqzJkzrj4zZsxQSkqK5syZo507dyoiIkJ9+/ZVcXGxx+cFAAB1h0dh5osvvtCtt96q5557ThEREXryySe1Y8eOKo+zdu1aDR8+XG3atFGHDh2UmpqqnJwc7dq1S9JPqzKzZ8/W5MmTNXjwYLVt21aLFi3S2bNntWTJknLHLCkpUVFRkdsGAADqLo/CTNu2bZWSkqKvv/5aqampys/P1x133KE2bdooJSVF3377rUfFFBYWSpIaNWok6afLWfn5+UpMTHT1cTqd6tmzp7Zt21buGMnJyQoLC3NtMTExHtUCAADscFU3APv7+2vQoEF6//339eqrr+rQoUMaP368oqOjNXToUOXl5VV6LGOMkpKSdMcdd6ht27aSpPz8fElSeHi4W9/w8HDXsV+aNGmSCgsLXVtubq6Hzw4AANjgqsJMRkaGRo8ercjISKWkpGj8+PE6dOiQNmzYoK+//lr33Xdfpcd66qmn9OWXX+q9994rc8zhcLjtG2PKtF3idDoVGhrqtgEAgLrL35MHpaSkKDU1VQcOHFD//v21ePFi9e/fX/Xq/ZSNmjdvrgULFuimm26q1HhPP/20PvzwQ23atEnR0dGu9oiICEk/rdBERka62gsKCsqs1gAAgH9PHq3MzJ8/X0OGDFFOTo5WrVqlX//6164gc8kNN9ygt99+u8JxjDF66qmntGLFCm3YsEHNmzd3O968eXNFREQoLS3N1Xb+/Hmlp6erW7dunpQOAADqGI9WZrKysq7YJzAwUMOGDauwz5gxY7RkyRL9/e9/V0hIiOs+mLCwMAUFBcnhcGjcuHGaNm2a4uLiFBcXp2nTpqlBgwYaMmSIJ6UDAIA6xqMwk5qaqoYNG+r+++93a//b3/6ms2fPXjHEXDJ//nxJUkJCQpnxhw8fLkl6/vnn9cMPP2j06NE6deqUOnfurPXr1yskJMST0gEAQB3j0WWm6dOnq0mTJmXamzZtqmnTplV6HGNMudulICP9dPPv1KlTlZeXp3Pnzik9Pd31bicAAACPwszRo0fL3N8iSbGxscrJybnqogAAACrLozDTtGlTffnll2Xav/jiCzVu3PiqiwIAAKgsj8LMgw8+qGeeeUaffvqpLl68qIsXL2rDhg0aO3asHnzwwequEQAA4LI8ugH4lVde0dGjR3XXXXfJ3/+nIUpLSzV06NAq3TMDAABwtTwKM4GBgVq2bJlefvllffHFFwoKClK7du0UGxtb3fUBAABUyKMwc0l8fLzi4+OrqxYAAIAq8yjMXLx4UQsXLtQnn3yigoIClZaWuh3fsGFDtRQHAABwJR6FmbFjx2rhwoW699571bZt28t+6SMAAEBN8yjMLF26VO+//7769+9f3fUAAABUiUdvzQ4MDFTLli2ruxYAAIAq8yjMPPfcc3rttddkjKnuegAAAKrEo8tMW7Zs0aeffqo1a9aoTZs2CggIcDu+YsWKaikOAADgSjwKM9dcc40GDRpU3bUAAABUmUdhJjU1tbrrAAAA8IhH98xI0oULF/Txxx9rwYIFKi4uliQdP35cp0+frrbiAAAArsSjlZmjR4/qnnvuUU5OjkpKStS3b1+FhIRoxowZOnfunN58883qrhMAAKBcHq3MjB07Vp06ddKpU6cUFBTkah80aJA++eSTaisOAADgSjx+N9PWrVsVGBjo1h4bG6uvv/66WgoDAACoDI9WZkpLS3Xx4sUy7ceOHVNISMhVFwUAAFBZHoWZvn37avbs2a59h8Oh06dPa8qUKXzFAQAAqFUeXWb64x//qF69eql169Y6d+6chgwZoqysLDVp0kTvvfdeddcIAABwWR6FmaioKGVmZuq9997T7t27VVpaqpEjR+rhhx92uyEYAACgpnkUZiQpKChII0aM0IgRI6qzHgAAgCrxKMwsXry4wuNDhw71qBgAAICq8ijMjB071m3/xx9/1NmzZxUYGKgGDRoQZgAAQK3x6N1Mp06dcttOnz6tAwcO6I477uAGYAAAUKs8/m6mX4qLi9P06dPLrNoAAADUpGoLM5Lk5+en48ePV+eQAAAAFfLonpkPP/zQbd8Yo7y8PM2ZM0fdu3evlsIAAAAqw6MwM3DgQLd9h8Oh6667Tr1799asWbOqoy4AAIBK8SjMlJaWVncdAAAAHqnWe2YAAABqm0crM0lJSZXum5KS4skpAAAAKsWjMLNnzx7t3r1bFy5cUKtWrSRJBw8elJ+fnzp27Ojq53A4qqdKAACAy/AozAwYMEAhISFatGiRrr32Wkk/fZDeY489pjvvvFPPPfdctRYJAABwOR7dMzNr1iwlJye7gowkXXvttXrllVd4NxMAAKhVHoWZoqIiffPNN2XaCwoKVFxcfNVFAQAAVJZHYWbQoEF67LHHtHz5ch07dkzHjh3T8uXLNXLkSA0ePLi6awQAALgsj+6ZefPNNzV+/Hg98sgj+vHHH38ayN9fI0eO1MyZM6u1QAAAgIp4FGYaNGigefPmaebMmTp06JCMMWrZsqWCg4Oruz4AAIAKXdWH5uXl5SkvL0/x8fEKDg6WMaa66gIAAKgUj8LMyZMndddddyk+Pl79+/dXXl6eJOnxxx/nbdkAAKBWeRRmnn32WQUEBCgnJ0cNGjRwtT/wwANau3ZtpcfZtGmTBgwYoKioKDkcDq1atcrt+PDhw+VwONy2Ll26eFIyAACoozy6Z2b9+vVat26doqOj3drj4uJ09OjRSo9z5swZdejQQY899pj+4z/+o9w+99xzj1JTU137gYGBnpQMAADqKI/CzJkzZ9xWZC45ceKEnE5npcfp16+f+vXrV2Efp9OpiIiIKtcIAAD+PXh0malHjx5avHixa9/hcKi0tFQzZ85Ur169qq04Sdq4caOaNm2q+Ph4PfHEEyooKKiwf0lJiYqKitw2AABQd3m0MjNz5kwlJCQoIyND58+f1/PPP6+9e/fqu+++09atW6utuH79+un+++9XbGyssrOz9T//8z/q3bu3du3addkVoOTkZL344ovVVgMAAPBtDuPh+6nz8/M1f/587dq1S6WlperYsaPGjBmjyMhIzwpxOLRy5UoNHDjwsn3y8vIUGxurpUuXXvaThktKSlRSUuLaLyoqUkxMjAoLCxUaGupRbRVpNvGjah9Tko5Mv7dGxgUAwAZFRUUKCwur1Ot3lVdmfvzxRyUmJmrBggW1vgISGRmp2NhYZWVlXbaP0+ms0n07AADAblW+ZyYgIED//Oc/5XA4aqKeCp08eVK5ubker/4AAIC6x6MbgIcOHaq33377qk9++vRpZWZmKjMzU5KUnZ2tzMxM5eTk6PTp0xo/frw+++wzHTlyRBs3btSAAQPUpEkTDRo06KrPDQAA6gaPbgA+f/68/vznPystLU2dOnUq851MKSkplRonIyPD7d1PSUlJkqRhw4Zp/vz5+uqrr7R48WJ9//33ioyMVK9evbRs2TKFhIR4UjYAAKiDqhRmDh8+rGbNmumf//ynOnbsKEk6ePCgW5+qXH5KSEio8Puc1q1bV5XyAADAv6EqhZm4uDjl5eXp008/lfTT1xe8/vrrCg8Pr5HiAAAArqRK98z8chVlzZo1OnPmTLUWBAAAUBUe3QB8iYcfUQMAAFBtqhRmLn1z9S/bAAAAvKVK98wYYzR8+HDXh9KdO3dOo0aNKvNuphUrVlRfhQAAABWoUpgZNmyY2/4jjzxSrcUAAABUVZXCTGpqak3VAQAA4JGrugEYAADA2wgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwmr+3C0DtazbxoxoZ98j0e2tkXAAAKsLKDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1bwaZjZt2qQBAwYoKipKDodDq1atcjtujNHUqVMVFRWloKAgJSQkaO/evd4pFgAA+CSvhpkzZ86oQ4cOmjNnTrnHZ8yYoZSUFM2ZM0c7d+5URESE+vbtq+Li4lquFAAA+Cqvfmhev3791K9fv3KPGWM0e/ZsTZ48WYMHD5YkLVq0SOHh4VqyZImefPLJ2iwVAAD4KJ+9ZyY7O1v5+flKTEx0tTmdTvXs2VPbtm277ONKSkpUVFTktgEAgLrLZ8NMfn6+JCk8PNytPTw83HWsPMnJyQoLC3NtMTExNVonAADwLp8NM5c4HA63fWNMmbafmzRpkgoLC11bbm5uTZcIAAC8yGe/aDIiIkLSTys0kZGRrvaCgoIyqzU/53Q65XQ6a7w+AADgG3x2ZaZ58+aKiIhQWlqaq+38+fNKT09Xt27dvFgZAADwJV5dmTl9+rT+7//+z7WfnZ2tzMxMNWrUSDfccIPGjRunadOmKS4uTnFxcZo2bZoaNGigIUOGeLFqAADgS7waZjIyMtSrVy/XflJSkiRp2LBhWrhwoZ5//nn98MMPGj16tE6dOqXOnTtr/fr1CgkJ8VbJAADAx3g1zCQkJMgYc9njDodDU6dO1dSpU2uvKAAAYBWfvWcGAACgMggzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKzm7+0CUL5mEz/ydgkAAFiBlRkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABW8/d2AUBd1WziRzUy7pHp99bIuABgK1ZmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACs5tNhZurUqXI4HG5bRESEt8sCAAA+xOc/Z6ZNmzb6+OOPXft+fn5erAYAAPganw8z/v7+VVqNKSkpUUlJiWu/qKioJsoCAAA+wqcvM0lSVlaWoqKi1Lx5cz344IM6fPhwhf2Tk5MVFhbm2mJiYmqpUgAA4A0+HWY6d+6sxYsXa926dXrrrbeUn5+vbt266eTJk5d9zKRJk1RYWOjacnNza7FiAABQ23z6MlO/fv1c/92uXTt17dpVN954oxYtWqSkpKRyH+N0OuV0OmurRAAA4GU+vTLzS8HBwWrXrp2ysrK8XQoAAPARVoWZkpIS7d+/X5GRkd4uBQAA+AifDjPjx49Xenq6srOz9fnnn+u3v/2tioqKNGzYMG+XBgAAfIRP3zNz7NgxPfTQQzpx4oSuu+46denSRdu3b1dsbKy3SwMAAD7Cp8PM0qVLvV0CAADwcT59mQkAAOBKCDMAAMBqhBkAAGA1n75nBnZpNvGjGhv7yPR7a2TcmqwZAFA7WJkBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKvxdQawAl87AAC4HFZmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDV+ARgAFaryU+HPjL93hobG0D1YWUGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKzG1xkAlqnJj++vKXwtQO2pqd8PfobwZazMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYzYowM2/ePDVv3lz169fXbbfdps2bN3u7JAAA4CN8PswsW7ZM48aN0+TJk7Vnzx7deeed6tevn3JycrxdGgAA8AE+H2ZSUlI0cuRIPf7447r55ps1e/ZsxcTEaP78+d4uDQAA+ACf/gTg8+fPa9euXZo4caJbe2JiorZt21buY0pKSlRSUuLaLywslCQVFRXVSI2lJWdrZFygLqmpvz+pZv8Ga7LumlJT82HjXMBul37njDFX7OvTYebEiRO6ePGiwsPD3drDw8OVn59f7mOSk5P14osvlmmPiYmpkRoBXFnYbG9X4Blb664JzAW8pbi4WGFhYRX28ekwc4nD4XDbN8aUabtk0qRJSkpKcu2Xlpbqu+++U+PGjS/7GE8VFRUpJiZGubm5Cg0Nrdax8f8xz7WDea4dzHPtYJ5rR03OszFGxcXFioqKumJfnw4zTZo0kZ+fX5lVmIKCgjKrNZc4nU45nU63tmuuuaamSpQkhYaG8sdSC5jn2sE81w7muXYwz7Wjpub5Sisyl/j0DcCBgYG67bbblJaW5taelpambt26eakqAADgS3x6ZUaSkpKS9Oijj6pTp07q2rWr/vSnPyknJ0ejRo3ydmkAAMAH+HyYeeCBB3Ty5Em99NJLysvLU9u2bbV69WrFxsZ6uzQ5nU5NmTKlzGUtVC/muXYwz7WDea4dzHPt8JV5dpjKvOcJAADAR/n0PTMAAABXQpgBAABWI8wAAACrEWYAAIDVCDMVmDdvnpo3b6769evrtttu0+bNmyvsn56erttuu03169dXixYt9Oabb9ZSpfarylyvWLFCffv21XXXXafQ0FB17dpV69atq8Vq7VXV3+lLtm7dKn9/f91yyy01W2AdUdV5Likp0eTJkxUbGyun06kbb7xR77zzTi1Va6+qzvO7776rDh06qEGDBoqMjNRjjz2mkydP1lK1dtq0aZMGDBigqKgoORwOrVq16oqP8cproUG5li5dagICAsxbb71l9u3bZ8aOHWuCg4PN0aNHy+1/+PBh06BBAzN27Fizb98+89Zbb5mAgACzfPnyWq7cPlWd67Fjx5pXX33V7Nixwxw8eNBMmjTJBAQEmN27d9dy5Xap6jxf8v3335sWLVqYxMRE06FDh9op1mKezPNvfvMb07lzZ5OWlmays7PN559/brZu3VqLVdunqvO8efNmU69ePfPaa6+Zw4cPm82bN5s2bdqYgQMH1nLldlm9erWZPHmy+eCDD4wks3Llygr7e+u1kDBzGbfffrsZNWqUW9tNN91kJk6cWG7/559/3tx0001ubU8++aTp0qVLjdVYV1R1rsvTunVr8+KLL1Z3aXWKp/P8wAMPmBdeeMFMmTKFMFMJVZ3nNWvWmLCwMHPy5MnaKK/OqOo8z5w507Ro0cKt7fXXXzfR0dE1VmNdU5kw463XQi4zleP8+fPatWuXEhMT3doTExO1bdu2ch/z2Weflel/9913KyMjQz/++GON1Wo7T+b6l0pLS1VcXKxGjRrVRIl1gqfznJqaqkOHDmnKlCk1XWKd4Mk8f/jhh+rUqZNmzJih66+/XvHx8Ro/frx++OGH2ijZSp7Mc7du3XTs2DGtXr1axhh98803Wr58ue69997aKPnfhrdeC33+E4C94cSJE7p48WKZL7MMDw8v86WXl+Tn55fb/8KFCzpx4oQiIyNrrF6beTLXvzRr1iydOXNGv/vd72qixDrBk3nOysrSxIkTtXnzZvn7809FZXgyz4cPH9aWLVtUv359rVy5UidOnNDo0aP13Xffcd/MZXgyz926ddO7776rBx54QOfOndOFCxf0m9/8Rm+88UZtlPxvw1uvhazMVMDhcLjtG2PKtF2pf3ntKKuqc33Je++9p6lTp2rZsmVq2rRpTZVXZ1R2ni9evKghQ4boxRdfVHx8fG2VV2dU5fe5tLRUDodD7777rm6//Xb1799fKSkpWrhwIaszV1CVed63b5+eeeYZ/eEPf9CuXbu0du1aZWdn8z1/NcAbr4X871Y5mjRpIj8/vzIJv6CgoEzivCQiIqLc/v7+/mrcuHGN1Wo7T+b6kmXLlmnkyJH629/+pj59+tRkmdar6jwXFxcrIyNDe/bs0VNPPSXppxddY4z8/f21fv169e7du1Zqt4knv8+RkZG6/vrrFRYW5mq7+eabZYzRsWPHFBcXV6M128iTeU5OTlb37t01YcIESVL79u0VHBysO++8U6+88gqr59XEW6+FrMyUIzAwULfddpvS0tLc2tPS0tStW7dyH9O1a9cy/devX69OnTopICCgxmq1nSdzLf20IjN8+HAtWbKEa96VUNV5Dg0N1VdffaXMzEzXNmrUKLVq1UqZmZnq3LlzbZVuFU9+n7t3767jx4/r9OnTrraDBw+qXr16io6OrtF6beXJPJ89e1b16rm/5Pn5+Un6/ysHuHpeey2s0duLLXbpbX9vv/222bdvnxk3bpwJDg42R44cMcYYM3HiRPPoo4+6+l96O9qzzz5r9u3bZ95++23eml1JVZ3rJUuWGH9/fzN37lyTl5fn2r7//ntvPQUrVHWef4l3M1VOVee5uLjYREdHm9/+9rdm7969Jj093cTFxZnHH3/cW0/BClWd59TUVOPv72/mzZtnDh06ZLZs2WI6depkbr/9dm89BSsUFxebPXv2mD179hhJJiUlxezZs8f1FnhfeS0kzFRg7ty5JjY21gQGBpqOHTua9PR017Fhw4aZnj17uvXfuHGjufXWW01gYKBp1qyZmT9/fi1XbK+qzHXPnj2NpDLbsGHDar9wy1T1d/rnCDOVV9V53r9/v+nTp48JCgoy0dHRJikpyZw9e7aWq7ZPVef59ddfN61btzZBQUEmMjLSPPzww+bYsWO1XLVdPv300wr/vfWV10KHMayvAQAAe3HPDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAK/YunWr2rVrp4CAAA0cOLDWz5+QkKBx48bV+nkBVD/CDGCB4cOHy+FwaPr06W7tq1atksPh8FJVVycpKUm33HKLsrOztXDhwlo//4oVK/Tyyy+79ps1a6bZs2fXeh0Arh5hBrBE/fr19eqrr+rUqVPeLqVaHDp0SL1791Z0dLSuueYaj8Y4f/68x+dv1KiRQkJCPH68L/jxxx+9XQLgEwgzgCX69OmjiIgIJScnV9jvgw8+UJs2beR0OtWsWTPNmjXL7XizZs00bdo0jRgxQiEhIbrhhhv0pz/9ya3PsWPH9OCDD6pRo0YKDg5Wp06d9Pnnn7uOz58/XzfeeKMCAwPVqlUr/eUvf3F7vMPh0J///GcNGjRIDRo0UFxcnD788ENJ0pEjR+RwOHTy5EmNGDFCDofDtTKTnp6u22+/XU6nU5GRkZo4caIuXLjgGjchIUFPPfWUkpKS1KRJE/Xt21cbN26Uw+HQunXrdOuttyooKEi9e/dWQUGB1qxZo5tvvlmhoaF66KGHdPbsWbexLl1mSkhI0NGjR/Xss8/K4XDI4XDozJkzCg0N1fLly92e2z/+8Q8FBweruLi43Plfvny52rVrp6CgIDVu3Fh9+vTRmTNnXMffeecd188nMjJSTz31lOtYTk6O7rvvPjVs2FChoaH63e9+p2+++cZ1fOrUqbrlllv0zjvvqEWLFnI6nTLGqLCwUP/5n/+ppk2bKjQ0VL1799YXX3xRbn1AnVTjX2UJ4KoNGzbM3HfffWbFihWmfv36Jjc31xhjzMqVK83P/4wzMjJMvXr1zEsvvWQOHDhgUlNTTVBQkElNTXX1iY2NNY0aNTJz5841WVlZJjk52dSrV8/s37/fGGNMcXGxadGihbnzzjvN5s2bTVZWllm2bJnZtm2bMcaYFStWmICAADN37lxz4MABM2vWLOPn52c2bNjgOockEx0dbZYsWWKysrLMM888Yxo2bGhOnjxpLly4YPLy8kxoaKiZPXu2ycvLM2fPnjXHjh0zDRo0MKNHjzb79+83K1euNE2aNDFTpkxxjduzZ0/TsGFDM2HCBPOvf/3L7N+/3/Wtvl26dDFbtmwxu3fvNi1btjQ9e/Y0iYmJZvfu3WbTpk2mcePGZvr06W5jjR071hhjzMmTJ010dLR56aWXTF5ensnLyzPGGPPEE0+Y/v37u/0sBg0aZIYOHVruz+n48ePG39/fpKSkmOzsbPPll1+auXPnmuLiYmOMMfPmzTP169c3s2fPNgcOHDA7duwwf/zjH40xxpSWlppbb73V3HHHHSYjI8Ns377ddOzY0e0biadMmWKCg4PN3XffbXbv3m2++OILU1paarp3724GDBhgdu7caQ4ePGiee+4507hxY3Py5Mkr/WoBdQJhBrDApTBjjDFdunQxI0aMMMaUDTNDhgwxffv2dXvshAkTTOvWrV37sbGx5pFHHnHtl5aWmqZNm5r58+cbY4xZsGCBCQkJuewLYbdu3cwTTzzh1nb//fe7vehLMi+88IJr//Tp08bhcJg1a9a42sLCwtxC1u9//3vTqlUrU1pa6mqbO3euadiwobl48aIx5qcAcsstt7id+1KY+fjjj11tycnJRpI5dOiQq+3JJ580d999t2v/52Hm0rxcChaXfP7558bPz898/fXXxhhjvv32WxMQEGA2btxY7tzs2rXLSDJHjhwp93hUVJSZPHlyucfWr19v/Pz8TE5Ojqtt7969RpLZsWOHMeanMBMQEGAKCgpcfT755BMTGhpqzp075zbejTfeaBYsWFDuuYC6hstMgGVeffVVLVq0SPv27StzbP/+/erevbtbW/fu3ZWVlaWLFy+62tq3b+/6b4fDoYiICBUUFEiSMjMzdeutt6pRo0blnv9y59i/f79b28/PERwcrJCQENc5Ljdu165d3W5o7t69u06fPq1jx4652jp16lTu439+vvDwcDVo0EAtWrRwa6vo/OW5/fbb1aZNGy1evFiS9Je//EU33HCDevToUW7/Dh066K677lK7du10//3366233nLd41RQUKDjx4/rrrvuKvex+/fvV0xMjGJiYlxtrVu31jXXXOM2t7Gxsbruuutc+7t27dLp06fVuHFjNWzY0LVlZ2fr0KFDVXq+gK0IM4BlevToobvvvlu///3vyxwzxpR5d5Mxpky/gIAAt32Hw6HS0lJJUlBQ0BVrKO8cv2yr6Bzlqaj2n7cHBweX+/ifn8/hcFT5/Jfz+OOPKzU1VZKUmpqqxx577LLvIPPz81NaWprWrFmj1q1b64033lCrVq2UnZ19xXkt7/mX1/7L519aWqrIyEhlZma6bQcOHNCECROq+nQBKxFmAAtNnz5d//jHP7Rt2za39tatW2vLli1ubdu2bVN8fLz8/PwqNXb79u2VmZmp7777rtzjN998c7nnuPnmm6vwDMpq3bq1tm3b5ha+tm3bppCQEF1//fVXNXZlBAYGuq1eXfLII48oJydHr7/+uvbu3athw4ZVOI7D4VD37t314osvas+ePQoMDNTKlSsVEhKiZs2a6ZNPPin3ca1bt1ZOTo5yc3Ndbfv27VNhYWGFc9uxY0fl5+fL399fLVu2dNuaNGlSyWcP2I0wA1ioXbt2evjhh/XGG2+4tT/33HP65JNP9PLLL+vgwYNatGiR5syZo/Hjx1d67IceekgREREaOHCgtm7dqsOHD+uDDz7QZ599JkmaMGGCFi5cqDfffFNZWVlKSUnRihUrqnSO8owePVq5ubl6+umn9a9//Ut///vfNWXKFCUlJalevZr/p6pZs2batGmTvv76a504ccLVfu2112rw4MGaMGGCEhMTFR0dfdkxPv/8c02bNk0ZGRnKycnRihUr9O2337rCyNSpUzVr1iy9/vrrysrK0u7du10/wz59+qh9+/Z6+OGHtXv3bu3YsUNDhw5Vz549L3tp7dLjunbtqoEDB2rdunU6cuSItm3bphdeeEEZGRnVNDuAbyPMAJZ6+eWXy1xC6tixo95//30tXbpUbdu21R/+8Ae99NJLGj58eKXHDQwM1Pr169W0aVP1799f7dq10/Tp010rOwMHDtRrr72mmTNnqk2bNlqwYIFSU1OVkJBwVc/n+uuv1+rVq7Vjxw516NBBo0aN0siRI/XCCy9c1biV9dJLL+nIkSO68cYb3e5JkaSRI0fq/PnzGjFiRIVjhIaGatOmTerfv7/i4+P1wgsvaNasWerXr58kadiwYZo9e7bmzZunNm3a6Ne//rWysrIk/bSis2rVKl177bXq0aOH+vTpoxYtWmjZsmUVntPhcGj16tXq0aOHRowYofj4eD344IM6cuSIwsPDr2JGAHs4THkX1AEALu+++67Gjh2r48ePKzAw0NvlAPgFf28XAAC+6uzZs8rOzlZycrKefPJJggzgo7jMBACXMWPGDN1yyy0KDw/XpEmTvF0OgMvgMhMAALAaKzMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNX+H8eC1ER7HxA2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.hist(nc_scores, bins=20, range=[0, 1])\n",
    "plt.xlabel(\"Nonconformity score\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "14521e74-11ea-40b6-b03a-b6fc837d65ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.972972972972973"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# apply finite sample correction\n",
    "\n",
    "n = y_calib.shape[0]\n",
    "\n",
    "q_level = np.ceil((n+1)*(confidence_level))/n\n",
    "q_level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "df26b84c-b07d-4a18-996f-703b486823c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5469768837244118"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# select the threshold such that 95% of the cases\n",
    "# are included (fall below)\n",
    "\n",
    "thres = np.quantile(nc_scores, q_level, method=\"higher\")\n",
    "thres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9adc520e-1afb-4971-895b-8cced5233bdc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzcElEQVR4nO3de3RU1f3+8WfIjSQkQUByMTHcEiTcRKkIqAQEFCgVaK2KlqtWCio0gF8otqC2BEVStApSqwFaEVoL1FauFQEBkbtaoJBCICCJEcRcJZFk//5gZX6OCSEZkszs+H6tddbi7LNnn8/sBOZhnzMzDmOMEQAAgKUaeLoAAACAq0GYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwmq+nC6htpaWlOnPmjEJCQuRwODxdDgAAqAJjjPLy8hQVFaUGDSpfe6n3YebMmTOKiYnxdBkAAMANp06dUnR0dKV96n2YCQkJkXRpMkJDQz1cDQBvVlBQoKioKEmX/iMUHBzs4YqA76/c3FzFxMQ4X8crU+/DTNmlpdDQUMIMgEr5+Pg4/xwaGkqYAbxAVW4R4QZgAABgNcIMAACwGmEGAABYrd7fMwMAqD9KS0tVXFzs6TJQA/z8/FzuU7sahBkAgBWKi4uVnp6u0tJST5eCGtK4cWNFRERc9efAEWYAAF7PGKPMzEz5+PgoJibmih+iBu9mjFFhYaGys7MlSZGRkVc1HmEGAOD1Ll68qMLCQkVFRSkoKMjT5aAGBAYGSpKys7PVvHnzq7rkRLQFAHi9kpISSZK/v7+HK0FNKgum33zzzVWNQ5gBAFiD79irX2rq50mYAQAAViPMAAAAqxFmAADwkBMnTsjhcOjAgQN1et7NmzfL4XDoq6++uqpxHA6HVq9efdnjdfX8CDMAANQCh8NR6TZq1ChPl1hv8NZsAABqQWZmpvPPK1as0G9+8xsdOXLE2RYYGKjz589Xe9ySkhI5HA4+a+dbmAkAgHWMMSooKPDIZoypUo0RERHOLSwsTA6Ho1xbmePHj6t3794KCgpS586d9eGHHzqPLV68WI0bN9a//vUvJSQkKCAgQCdPnlRxcbGefPJJXXfddQoODla3bt20efNm5+NOnjypwYMH65prrlFwcLDat2+vNWvWuNS4d+9ede3aVUFBQerRo4dL2JKkhQsXqnXr1vL391fbtm315z//udLnvGvXLnXp0kUNGzZU165dtX///irN1dViZeYqtZj2bq2Me2LOoFoZFwDqg8LCQjVq1Mgj587Pz1dwcHCNjjljxgy98MILiouL04wZM/TAAw/of//7n3x9L71MFxYWKjk5WX/605/UtGlTNW/eXKNHj9aJEye0fPlyRUVFadWqVbr77rv16aefKi4uThMmTFBxcbG2bt2q4OBgHTp0qNyczZgxQ/PmzdO1116rcePGacyYMdq+fbskadWqVZo4caLmz5+vvn376l//+pdGjx6t6Oho9e7du9xzKCgo0A9/+EP16dNHf/nLX5Senq6JEyfW6DxdDmEGAAAPmzJligYNuvSf2Kefflrt27fX//73P91www2SLn2o3IIFC9S5c2dJ0rFjx/TWW2/p9OnTioqKco6xbt06paamavbs2crIyNCPf/xjdezYUZLUqlWrcuf93e9+p169ekmSpk2bpkGDBunChQtq2LChXnjhBY0aNUrjx4+XJCUlJWnnzp164YUXKgwzb775pkpKSvTGG28oKChI7du31+nTp/WLX/yihmerPMIMAMA6QUFBys/P99i5a1qnTp2cfy77nqLs7GxnmPH393fps2/fPhljFB8f7zJOUVGRmjZtKkl64okn9Itf/EIbNmxQ37599eMf/9hljMrOe/311+vw4cP6+c9/7tK/Z8+eevHFFyt8DocPH1bnzp1d5qd79+5Vm4CrRJgBAFjH4XDU+KUeT/Lz83P+uexTcb/97eCBgYEun5ZbWloqHx8f7d27t9x3GpVdSnr44Yd111136d1339WGDRuUnJysefPm6fHHH6/yeb/7Cb3GmMt+am9V7yWqDdwADACAZbp06aKSkhJlZ2erTZs2LltERISzX0xMjMaNG6eVK1dq8uTJeu2116p8jnbt2mnbtm0ubTt27FC7du0q7J+QkKCPP/5YX3/9tbNt586d1Xxm7iHMAABgmfj4eD344IMaMWKEVq5cqfT0dO3evVvPPfec8x1LkyZN0vr165Wenq59+/Zp06ZNlw0iFZk6daoWL16sV199VWlpaUpJSdHKlSs1ZcqUCvsPHz5cDRo00NixY3Xo0CGtWbNGL7zwQo083yshzAAAYKHU1FSNGDFCkydPVtu2bfWjH/1IH330kWJiYiRd+jyaCRMmqF27drr77rvVtm1bLViwoMrjDxkyRC+++KLmzp2r9u3ba9GiRUpNTVViYmKF/Rs1aqR//vOfOnTokLp06aIZM2boueeeq4mnekUO48mLXHUgNzdXYWFhysnJUWhoaI2Pz1uzgfqjoKDAeb9Bbbz9Fu67cOGC0tPT1bJlSzVs2NDT5aCGVPZzrc7rNyszAADAaoQZAABgNcIMAACwGmEGAGCNen6b5/dOTf08CTMAAK9X9sFwxcXFHq4ENamwsFCS64f3uYNPAAYAeD1fX18FBQXpiy++kJ+fnxo04P/iNjPGqLCwUNnZ2WrcuHG5TzGuLsIMAMDrORwORUZGKj09XSdPnvR0OaghjRs3dvnEYncRZgAAVvD391dcXByXmuoJPz+/q16RKUOYAQBYo0GDBnxoHsrhoiMAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArObRMLNw4UJ16tRJoaGhCg0NVffu3bV27VrncWOMZs2apaioKAUGBioxMVEHDx70YMUAAMDbeDTMREdHa86cOdqzZ4/27NmjPn366J577nEGlueff14pKSl6+eWXtXv3bkVERKhfv37Ky8vzZNkAAMCLeDTMDB48WAMHDlR8fLzi4+P1u9/9To0aNdLOnTtljNH8+fM1Y8YMDRs2TB06dNCSJUtUWFioZcuWebJsAADgRbzmnpmSkhItX75cBQUF6t69u9LT05WVlaX+/fs7+wQEBKhXr17asWPHZccpKipSbm6uywYAAOovj4eZTz/9VI0aNVJAQIDGjRunVatWKSEhQVlZWZKk8PBwl/7h4eHOYxVJTk5WWFiYc4uJianV+gEAgGd5PMy0bdtWBw4c0M6dO/WLX/xCI0eO1KFDh5zHHQ6HS39jTLm2b5s+fbpycnKc26lTp2qtdgAA4Hm+ni7A399fbdq0kSR17dpVu3fv1osvvqj/+7//kyRlZWUpMjLS2T87O7vcas23BQQEKCAgoHaLBgAAXsPjKzPfZYxRUVGRWrZsqYiICG3cuNF5rLi4WFu2bFGPHj08WCEAAPAmHl2Z+dWvfqUBAwYoJiZGeXl5Wr58uTZv3qx169bJ4XBo0qRJmj17tuLi4hQXF6fZs2crKChIw4cP92TZAADAi3g0zHz++ef62c9+pszMTIWFhalTp05at26d+vXrJ0l68skn9fXXX2v8+PE6f/68unXrpg0bNigkJMSTZQMAAC/iMMYYTxdRm3JzcxUWFqacnByFhobW+Pgtpr1b42NK0ok5g2plXACXV1BQoEaNGkmS8vPzFRwc7OGKgO+v6rx+e909MwAAANVBmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABW82iYSU5O1g9+8AOFhISoefPmGjJkiI4cOeLSZ9SoUXI4HC7brbfe6qGKAQCAt/FomNmyZYsmTJignTt3auPGjbp48aL69++vgoICl3533323MjMznduaNWs8VDEAAPA2vp48+bp161z2U1NT1bx5c+3du1d33HGHsz0gIEARERF1XR4AALCAV90zk5OTI0lq0qSJS/vmzZvVvHlzxcfH65FHHlF2dvZlxygqKlJubq7LBgAA6i+vCTPGGCUlJem2225Thw4dnO0DBgzQm2++qU2bNmnevHnavXu3+vTpo6KiogrHSU5OVlhYmHOLiYmpq6cAAAA8wGGMMZ4uQpImTJigd999V9u2bVN0dPRl+2VmZio2NlbLly/XsGHDyh0vKipyCTq5ubmKiYlRTk6OQkNDa7zuFtPerfExJenEnEG1Mi6AyysoKFCjRo0kSfn5+QoODvZwRcD3V25ursLCwqr0+u3Re2bKPP7443rnnXe0devWSoOMJEVGRio2NlZpaWkVHg8ICFBAQEBtlAkAALyQR8OMMUaPP/64Vq1apc2bN6tly5ZXfMy5c+d06tQpRUZG1kGFAADA23n0npkJEyboL3/5i5YtW6aQkBBlZWUpKytLX3/9taRLy7xTpkzRhx9+qBMnTmjz5s0aPHiwmjVrpqFDh3qydAAA4CU8ujKzcOFCSVJiYqJLe2pqqkaNGiUfHx99+umnWrp0qb766itFRkaqd+/eWrFihUJCQjxQMQAA8DYev8xUmcDAQK1fv76OqgEAADbymrdmAwAAuIMwAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFZzK8ykp6fXyMmTk5P1gx/8QCEhIWrevLmGDBmiI0eOuPQxxmjWrFmKiopSYGCgEhMTdfDgwRo5PwAAsJ9bYaZNmzbq3bu3/vKXv+jChQtun3zLli2aMGGCdu7cqY0bN+rixYvq37+/CgoKnH2ef/55paSk6OWXX9bu3bsVERGhfv36KS8vz+3zAgCA+sOtMPPxxx+rS5cumjx5siIiIvToo49q165d1R5n3bp1GjVqlNq3b6/OnTsrNTVVGRkZ2rt3r6RLqzLz58/XjBkzNGzYMHXo0EFLlixRYWGhli1bVuGYRUVFys3NddkAAED95VaY6dChg1JSUvTZZ58pNTVVWVlZuu2229S+fXulpKToiy++cKuYnJwcSVKTJk0kXbqclZWVpf79+zv7BAQEqFevXtqxY0eFYyQnJyssLMy5xcTEuFULAACww1XdAOzr66uhQ4fqr3/9q5577jkdO3ZMU6ZMUXR0tEaMGKHMzMwqj2WMUVJSkm677TZ16NBBkpSVlSVJCg8Pd+kbHh7uPPZd06dPV05OjnM7deqUm88OAADY4KrCzJ49ezR+/HhFRkYqJSVFU6ZM0bFjx7Rp0yZ99tlnuueee6o81mOPPaZPPvlEb731VrljDofDZd8YU66tTEBAgEJDQ102AABQf/m686CUlBSlpqbqyJEjGjhwoJYuXaqBAweqQYNL2ahly5ZatGiRbrjhhiqN9/jjj+udd97R1q1bFR0d7WyPiIiQdGmFJjIy0tmenZ1dbrUGAAB8P7m1MrNw4UINHz5cGRkZWr16tX74wx86g0yZ66+/Xq+//nql4xhj9Nhjj2nlypXatGmTWrZs6XK8ZcuWioiI0MaNG51txcXF2rJli3r06OFO6QAAoJ5xa2UmLS3tin38/f01cuTISvtMmDBBy5Yt0z/+8Q+FhIQ474MJCwtTYGCgHA6HJk2apNmzZysuLk5xcXGaPXu2goKCNHz4cHdKBwAA9YxbYSY1NVWNGjXSvffe69L+t7/9TYWFhVcMMWUWLlwoSUpMTCw3/qhRoyRJTz75pL7++muNHz9e58+fV7du3bRhwwaFhIS4UzoAAKhn3LrMNGfOHDVr1qxce/PmzTV79uwqj2OMqXArCzLSpZt/Z82apczMTF24cEFbtmxxvtsJAADArTBz8uTJcve3SFJsbKwyMjKuuigAAICqcivMNG/eXJ988km59o8//lhNmza96qIAAACqyq0wc//99+uJJ57Q+++/r5KSEpWUlGjTpk2aOHGi7r///pquEQAA4LLcugH4t7/9rU6ePKk777xTvr6XhigtLdWIESOqdc8MAADA1XIrzPj7+2vFihV69tln9fHHHyswMFAdO3ZUbGxsTdcHAABQKbfCTJn4+HjFx8fXVC0AAADV5laYKSkp0eLFi/Xee+8pOztbpaWlLsc3bdpUI8UBAABciVthZuLEiVq8eLEGDRqkDh06XPZLHwEAAGqbW2Fm+fLl+utf/6qBAwfWdD0AAADV4tZbs/39/dWmTZuargUAAKDa3AozkydP1osvvihjTE3XAwAAUC1uXWbatm2b3n//fa1du1bt27eXn5+fy/GVK1fWSHEAAABX4laYady4sYYOHVrTtQAAAFSbW2EmNTW1pusAAABwi1v3zEjSxYsX9e9//1uLFi1SXl6eJOnMmTPKz8+vseIAAACuxK2VmZMnT+ruu+9WRkaGioqK1K9fP4WEhOj555/XhQsX9Oqrr9Z0nQAAABVya2Vm4sSJ6tq1q86fP6/AwEBn+9ChQ/Xee+/VWHEAAABX4va7mbZv3y5/f3+X9tjYWH322Wc1UhgAAEBVuLUyU1paqpKSknLtp0+fVkhIyFUXBQAAUFVuhZl+/fpp/vz5zn2Hw6H8/HzNnDmTrzgAAAB1yq3LTL///e/Vu3dvJSQk6MKFCxo+fLjS0tLUrFkzvfXWWzVdIwAAwGW5FWaioqJ04MABvfXWW9q3b59KS0s1duxYPfjggy43BAMAANQ2t8KMJAUGBmrMmDEaM2ZMTdYDAABQLW6FmaVLl1Z6fMSIEW4VAwAAUF1uhZmJEye67H/zzTcqLCyUv7+/goKCCDMAAKDOuPVupvPnz7ts+fn5OnLkiG677TZuAAYAAHXK7e9m+q64uDjNmTOn3KoNAABAbaqxMCNJPj4+OnPmTE0OCQAAUCm37pl55513XPaNMcrMzNTLL7+snj171khhAAAAVeFWmBkyZIjLvsPh0LXXXqs+ffpo3rx5NVEXAABAlbgVZkpLS2u6DgAAALfU6D0zAAAAdc2tlZmkpKQq901JSXHnFAAAAFXiVpjZv3+/9u3bp4sXL6pt27aSpKNHj8rHx0c33XSTs5/D4aiZKgEAAC7DrTAzePBghYSEaMmSJbrmmmskXfogvdGjR+v222/X5MmTa7RIAACAy3Hrnpl58+YpOTnZGWQk6ZprrtFvf/tb3s0EAADqlFthJjc3V59//nm59uzsbOXl5V11UQAAAFXlVpgZOnSoRo8erbffflunT5/W6dOn9fbbb2vs2LEaNmxYTdcIAABwWW7dM/Pqq69qypQpeuihh/TNN99cGsjXV2PHjtXcuXNrtEAAAIDKuBVmgoKCtGDBAs2dO1fHjh2TMUZt2rRRcHBwTdcHAABQqav60LzMzExlZmYqPj5ewcHBMsbUVF0AAABV4laYOXfunO68807Fx8dr4MCByszMlCQ9/PDDvC0bAADUKbfCzC9/+Uv5+fkpIyNDQUFBzvb77rtP69atq/I4W7du1eDBgxUVFSWHw6HVq1e7HB81apQcDofLduutt7pTMgAAqKfcumdmw4YNWr9+vaKjo13a4+LidPLkySqPU1BQoM6dO2v06NH68Y9/XGGfu+++W6mpqc59f39/d0oGAAD1lFthpqCgwGVFpszZs2cVEBBQ5XEGDBigAQMGVNonICBAERER1a4RAAB8P7h1memOO+7Q0qVLnfsOh0OlpaWaO3euevfuXWPFSdLmzZvVvHlzxcfH65FHHlF2dnal/YuKipSbm+uyAQCA+sutlZm5c+cqMTFRe/bsUXFxsZ588kkdPHhQX375pbZv315jxQ0YMED33nuvYmNjlZ6erl//+tfq06eP9u7de9kVoOTkZD399NM1VgMAAPBubq3MJCQk6JNPPtEtt9yifv36qaCgQMOGDdP+/fvVunXrGivuvvvu06BBg9ShQwcNHjxYa9eu1dGjR/Xuu+9e9jHTp09XTk6Oczt16lSN1QMAALxPtVdmvvnmG/Xv31+LFi2q8xWQyMhIxcbGKi0t7bJ9AgICqnXfDgAAsFu1V2b8/Pz0n//8Rw6HozbqqdS5c+d06tQpRUZG1vm5AQCAd3LrMtOIESP0+uuvX/XJ8/PzdeDAAR04cECSlJ6ergMHDigjI0P5+fmaMmWKPvzwQ504cUKbN2/W4MGD1axZMw0dOvSqzw0AAOoHt24ALi4u1p/+9Cdt3LhRXbt2LfedTCkpKVUaZ8+ePS7vfkpKSpIkjRw5UgsXLtSnn36qpUuX6quvvlJkZKR69+6tFStWKCQkxJ2yAQBAPVStMHP8+HG1aNFC//nPf3TTTTdJko4ePerSpzqXnxITEyv9Pqf169dXpzwAAPA9VK0wExcXp8zMTL3//vuSLr3b6KWXXlJ4eHitFAcAAHAl1bpn5rurKGvXrlVBQUGNFgQAAFAdbt0AXKayS0QAAAB1oVphpuybq7/bBgAA4CnVumfGGKNRo0Y5P5TuwoULGjduXLl3M61cubLmKgQAAKhEtcLMyJEjXfYfeuihGi0GAACguqoVZlJTU2urDgAAALdc1Q3AAAAAnkaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1Xw9XQDqXotp79bKuCfmDKqVcQEAqAwrMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAVvNomNm6dasGDx6sqKgoORwOrV692uW4MUazZs1SVFSUAgMDlZiYqIMHD3qmWAAA4JU8GmYKCgrUuXNnvfzyyxUef/7555WSkqKXX35Zu3fvVkREhPr166e8vLw6rhQAAHgrj35o3oABAzRgwIAKjxljNH/+fM2YMUPDhg2TJC1ZskTh4eFatmyZHn300bosFQAAeCmvvWcmPT1dWVlZ6t+/v7MtICBAvXr10o4dOy77uKKiIuXm5rpsAACg/vLaMJOVlSVJCg8Pd2kPDw93HqtIcnKywsLCnFtMTEyt1gkAADzLa8NMGYfD4bJvjCnX9m3Tp09XTk6Oczt16lRtlwgAADzIa79oMiIiQtKlFZrIyEhne3Z2drnVmm8LCAhQQEBArdcHAAC8g9euzLRs2VIRERHauHGjs624uFhbtmxRjx49PFgZAADwJh5dmcnPz9f//vc/5356eroOHDigJk2a6Prrr9ekSZM0e/ZsxcXFKS4uTrNnz1ZQUJCGDx/uwaoBAIA38WiY2bNnj3r37u3cT0pKkiSNHDlSixcv1pNPPqmvv/5a48eP1/nz59WtWzdt2LBBISEhnioZAAB4GY+GmcTERBljLnvc4XBo1qxZmjVrVt0VBQAArOK198wAAABUBWEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFiNMAMAAKxGmAEAAFYjzAAAAKsRZgAAgNV8PV0AKtZi2rueLgEAACuwMgMAAKxGmAEAAFYjzAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqvp4uAKivWkx7t1bGPTFnUK2MCwC2YmUGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDVCDMAAMBqXh1mZs2aJYfD4bJFRER4uiwAAOBFvP5zZtq3b69///vfzn0fHx8PVgMAALyN14cZX1/faq3GFBUVqaioyLmfm5tbG2UBAAAv4dWXmSQpLS1NUVFRatmype6//34dP3680v7JyckKCwtzbjExMXVUKQAA8ASvDjPdunXT0qVLtX79er322mvKyspSjx49dO7cucs+Zvr06crJyXFup06dqsOKAQBAXfPqy0wDBgxw/rljx47q3r27WrdurSVLligpKanCxwQEBCggIKCuSgQAAB7m1Ssz3xUcHKyOHTsqLS3N06UAAAAvYVWYKSoq0uHDhxUZGenpUgAAgJfw6jAzZcoUbdmyRenp6froo4/0k5/8RLm5uRo5cqSnSwMAAF7Cq++ZOX36tB544AGdPXtW1157rW699Vbt3LlTsbGxni4NAAB4Ca8OM8uXL/d0CQAAwMt59WUmAACAKyHMAAAAqxFmAACA1bz6nhnYpcW0d2tt7BNzBtXKuLVZMwCgbrAyAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABW4+sMYAW+dgAAcDmszAAAAKsRZgAAgNUIMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAq/EJwACsVpOfDl1afMH553a/XqeMlB/X2NgAag8rMwAAwGqEGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAABgNb7OALBMTX58f105MWeQp0v43qit3w9+hvBmrMwAAACrEWYAAIDVCDMAAMBqhBkAAGA1wgwAALAaYQYAAFjNijCzYMECtWzZUg0bNtTNN9+sDz74wNMlAQAAL+H1YWbFihWaNGmSZsyYof379+v222/XgAEDlJGR4enSAACAF/D6MJOSkqKxY8fq4YcfVrt27TR//nzFxMRo4cKFni4NAAB4Aa/+BODi4mLt3btX06ZNc2nv37+/duzYUeFjioqKVFRU5NzPycmRJOXm5tZKjaVFhbUyLlCf1NbfP6lm/w6WFl9wGbc2664ttfVvko1zAbuV/c4ZY67Y16vDzNmzZ1VSUqLw8HCX9vDwcGVlZVX4mOTkZD399NPl2mNiYmqlRgBXFjbf0xVU32cLRihsgaer8B42/gxRP+Tl5SksLKzSPl4dZso4HA6XfWNMubYy06dPV1JSknO/tLRUX375pZo2bXrZx7grNzdXMTExOnXqlEJDQ2t0bPx/zHPdYJ7rBvNcN5jnulGb82yMUV5enqKioq7Y16vDTLNmzeTj41NuFSY7O7vcak2ZgIAABQQEuLQ1bty4tkqUJIWGhvKXpQ4wz3WDea4bzHPdYJ7rRm3N85VWZMp49Q3A/v7+uvnmm7Vx40aX9o0bN6pHjx4eqgoAAHgTr16ZkaSkpCT97Gc/U9euXdW9e3f98Y9/VEZGhsaNG+fp0gAAgBfw+jBz33336dy5c3rmmWeUmZmpDh06aM2aNYqNjfV0aQoICNDMmTPLXdZCzWKe6wbzXDeY57rBPNcNb5lnh6nKe54AAAC8lFffMwMAAHAlhBkAAGA1wgwAALAaYQYAAFiNMFOJBQsWqGXLlmrYsKFuvvlmffDBB5X237Jli26++WY1bNhQrVq10quvvlpHldqvOnO9cuVK9evXT9dee61CQ0PVvXt3rV+/vg6rtVd1f6fLbN++Xb6+vrrxxhtrt8B6orrzXFRUpBkzZig2NlYBAQFq3bq13njjjTqq1l7Vnec333xTnTt3VlBQkCIjIzV69GidO3eujqq109atWzV48GBFRUXJ4XBo9erVV3yMR14LDSq0fPly4+fnZ1577TVz6NAhM3HiRBMcHGxOnjxZYf/jx4+boKAgM3HiRHPo0CHz2muvGT8/P/P222/XceX2qe5cT5w40Tz33HNm165d5ujRo2b69OnGz8/P7Nu3r44rt0t157nMV199ZVq1amX69+9vOnfuXDfFWsydef7Rj35kunXrZjZu3GjS09PNRx99ZLZv316HVdunuvP8wQcfmAYNGpgXX3zRHD9+3HzwwQemffv2ZsiQIXVcuV3WrFljZsyYYf7+978bSWbVqlWV9vfUayFh5jJuueUWM27cOJe2G264wUybNq3C/k8++aS54YYbXNoeffRRc+utt9ZajfVFdee6IgkJCebpp5+u6dLqFXfn+b777jNPPfWUmTlzJmGmCqo7z2vXrjVhYWHm3LlzdVFevVHdeZ47d65p1aqVS9tLL71koqOja63G+qYqYcZTr4VcZqpAcXGx9u7dq/79+7u09+/fXzt27KjwMR9++GG5/nfddZf27Nmjb775ptZqtZ07c/1dpaWlysvLU5MmTWqjxHrB3XlOTU3VsWPHNHPmzNousV5wZ57feecdde3aVc8//7yuu+46xcfHa8qUKfr666/romQruTPPPXr00OnTp7VmzRoZY/T555/r7bff1qBBg+qi5O8NT70Wev0nAHvC2bNnVVJSUu7LLMPDw8t96WWZrKysCvtfvHhRZ8+eVWRkZK3VazN35vq75s2bp4KCAv30pz+tjRLrBXfmOS0tTdOmTdMHH3wgX1/+qagKd+b5+PHj2rZtmxo2bKhVq1bp7NmzGj9+vL788kvum7kMd+a5R48eevPNN3XffffpwoULunjxon70ox/pD3/4Q12U/L3hqddCVmYq4XA4XPaNMeXartS/onaUV925LvPWW29p1qxZWrFihZo3b15b5dUbVZ3nkpISDR8+XE8//bTi4+Prqrx6ozq/z6WlpXI4HHrzzTd1yy23aODAgUpJSdHixYtZnbmC6szzoUOH9MQTT+g3v/mN9u7dq3Xr1ik9PZ3v+asFnngt5L9bFWjWrJl8fHzKJfzs7OxyibNMREREhf19fX3VtGnTWqvVdu7MdZkVK1Zo7Nix+tvf/qa+ffvWZpnWq+485+Xlac+ePdq/f78ee+wxSZdedI0x8vX11YYNG9SnT586qd0m7vw+R0ZG6rrrrlNYWJizrV27djLG6PTp04qLi6vVmm3kzjwnJyerZ8+emjp1qiSpU6dOCg4O1u23367f/va3rJ7XEE+9FrIyUwF/f3/dfPPN2rhxo0v7xo0b1aNHjwof071793L9N2zYoK5du8rPz6/WarWdO3MtXVqRGTVqlJYtW8Y17yqo7jyHhobq008/1YEDB5zbuHHj1LZtWx04cEDdunWrq9Kt4s7vc8+ePXXmzBnl5+c7244ePaoGDRooOjq6Vuu1lTvzXFhYqAYNXF/yfHx8JP3/lQNcPY+9Ftbq7cUWK3vb3+uvv24OHTpkJk2aZIKDg82JEyeMMcZMmzbN/OxnP3P2L3s72i9/+Utz6NAh8/rrr/PW7Cqq7lwvW7bM+Pr6mldeecVkZmY6t6+++spTT8EK1Z3n7+LdTFVT3XnOy8sz0dHR5ic/+Yk5ePCg2bJli4mLizMPP/ywp56CFao7z6mpqcbX19csWLDAHDt2zGzbts107drV3HLLLZ56ClbIy8sz+/fvN/v37zeSTEpKitm/f7/zLfDe8lpImKnEK6+8YmJjY42/v7+56aabzJYtW5zHRo4caXr16uXSf/PmzaZLly7G39/ftGjRwixcuLCOK7ZXdea6V69eRlK5beTIkXVfuGWq+zv9bYSZqqvuPB8+fNj07dvXBAYGmujoaJOUlGQKCwvruGr7VHeeX3rpJZOQkGACAwNNZGSkefDBB83p06fruGq7vP/++5X+e+str4UOY1hfAwAA9uKeGQAAYDXCDAAAsBphBgAAWI0wAwAArEaYAQAAViPMAAAAqxFmAACA1QgzAADAaoQZAB6xfft2dezYUX5+fhoyZEidnz8xMVGTJk2q8/MCqHmEGcACo0aNksPh0Jw5c1zaV69eLYfD4aGqrk5SUpJuvPFGpaena/HixXV+/pUrV+rZZ5917rdo0ULz58+v8zoAXD3CDGCJhg0b6rnnntP58+c9XUqNOHbsmPr06aPo6Gg1btzYrTGKi4vdPn+TJk0UEhLi9uO9wTfffOPpEgCvQJgBLNG3b19FREQoOTm50n5///vf1b59ewUEBKhFixaaN2+ey/EWLVpo9uzZGjNmjEJCQnT99dfrj3/8o0uf06dP6/7771eTJk0UHBysrl276qOPPnIeX7hwoVq3bi1/f3+1bdtWf/7zn10e73A49Kc//UlDhw5VUFCQ4uLi9M4770iSTpw4IYfDoXPnzmnMmDFyOBzOlZktW7bolltuUUBAgCIjIzVt2jRdvHjROW5iYqIee+wxJSUlqVmzZurXr582b94sh8Oh9evXq0uXLgoMDFSfPn2UnZ2ttWvXql27dgoNDdUDDzygwsJCl7HKLjMlJibq5MmT+uUvfymHwyGHw6GCggKFhobq7bffdnlu//znPxUcHKy8vLwK5//tt99Wx44dFRgYqKZNm6pv374qKChwHn/jjTecP5/IyEg99thjzmMZGRm655571KhRI4WGhuqnP/2pPv/8c+fxWbNm6cYbb9Qbb7yhVq1aKSAgQMYY5eTk6Oc//7maN2+u0NBQ9enTRx9//HGF9QH1Uq1/lSWAqzZy5Ehzzz33mJUrV5qGDRuaU6dOGWOMWbVqlfn2X+M9e/aYBg0amGeeecYcOXLEpKammsDAQJOamursExsba5o0aWJeeeUVk5aWZpKTk02DBg3M4cOHjTHG5OXlmVatWpnbb7/dfPDBByYtLc2sWLHC7NixwxhjzMqVK42fn5955ZVXzJEjR8y8efOMj4+P2bRpk/Mckkx0dLRZtmyZSUtLM0888YRp1KiROXfunLl48aLJzMw0oaGhZv78+SYzM9MUFhaa06dPm6CgIDN+/Hhz+PBhs2rVKtOsWTMzc+ZM57i9evUyjRo1MlOnTjX//e9/zeHDh53f6nvrrbeabdu2mX379pk2bdqYXr16mf79+5t9+/aZrVu3mqZNm5o5c+a4jDVx4kRjjDHnzp0z0dHR5plnnjGZmZkmMzPTGGPMI488YgYOHOjysxg6dKgZMWJEhT+nM2fOGF9fX5OSkmLS09PNJ598Yl555RWTl5dnjDFmwYIFpmHDhmb+/PnmyJEjZteuXeb3v/+9McaY0tJS06VLF3PbbbeZPXv2mJ07d5qbbrrJ5RuJZ86caYKDg81dd91l9u3bZz7++GNTWlpqevbsaQYPHmx2795tjh49aiZPnmyaNm1qzp07d6VfLaBeIMwAFigLM8YYc+utt5oxY8YYY8qHmeHDh5t+/fq5PHbq1KkmISHBuR8bG2seeugh535paalp3ry5WbhwoTHGmEWLFpmQkJDLvhD26NHDPPLIIy5t9957r8uLviTz1FNPOffz8/ONw+Ewa9eudbaFhYW5hKxf/epXpm3btqa0tNTZ9sorr5hGjRqZkpISY8ylAHLjjTe6nLsszPz73/92tiUnJxtJ5tixY862Rx991Nx1113O/W+HmbJ5KQsWZT766CPj4+NjPvvsM2OMMV988YXx8/MzmzdvrnBu9u7daySZEydOVHg8KirKzJgxo8JjGzZsMD4+PiYjI8PZdvDgQSPJ7Nq1yxhzKcz4+fmZ7OxsZ5/33nvPhIaGmgsXLriM17p1a7No0aIKzwXUN1xmAizz3HPPacmSJTp06FC5Y4cPH1bPnj1d2nr27Km0tDSVlJQ42zp16uT8s8PhUEREhLKzsyVJBw4cUJcuXdSkSZMKz3+5cxw+fNil7dvnCA4OVkhIiPMclxu3e/fuLjc09+zZU/n5+Tp9+rSzrWvXrhU+/tvnCw8PV1BQkFq1auXSVtn5K3LLLbeoffv2Wrp0qSTpz3/+s66//nrdcccdFfbv3Lmz7rzzTnXs2FH33nuvXnvtNec9TtnZ2Tpz5ozuvPPOCh97+PBhxcTEKCYmxtmWkJCgxo0bu8xtbGysrr32Wuf+3r17lZ+fr6ZNm6pRo0bOLT09XceOHavW8wVsRZgBLHPHHXforrvu0q9+9atyx4wx5d7dZIwp18/Pz89l3+FwqLS0VJIUGBh4xRoqOsd32yo7R0Uqq/3b7cHBwRU+/tvnczgc1T7/5Tz88MNKTU2VJKWmpmr06NGXfQeZj4+PNm7cqLVr1yohIUF/+MMf1LZtW6Wnp19xXit6/hW1f/f5l5aWKjIyUgcOHHDZjhw5oqlTp1b36QJWIswAFpozZ47++c9/aseOHS7tCQkJ2rZtm0vbjh07FB8fLx8fnyqN3alTJx04cEBffvllhcfbtWtX4TnatWtXjWdQXkJCgnbs2OESvnbs2KGQkBBdd911VzV2Vfj7+7usXpV56KGHlJGRoZdeekkHDx7UyJEjKx3H4XCoZ8+eevrpp7V//375+/tr1apVCgkJUYsWLfTee+9V+LiEhARlZGTo1KlTzrZDhw4pJyen0rm96aablJWVJV9fX7Vp08Zla9asWRWfPWA3wgxgoY4dO+rBBx/UH/7wB5f2yZMn67333tOzzz6ro0ePasmSJXr55Zc1ZcqUKo/9wAMPKCIiQkOGDNH27dt1/Phx/f3vf9eHH34oSZo6daoWL16sV199VWlpaUpJSdHKlSurdY6KjB8/XqdOndLjjz+u//73v/rHP/6hmTNnKikpSQ0a1P4/VS1atNDWrVv12Wef6ezZs872a665RsOGDdPUqVPVv39/RUdHX3aMjz76SLNnz9aePXuUkZGhlStX6osvvnCGkVmzZmnevHl66aWXlJaWpn379jl/hn379lWnTp304IMPat++fdq1a5dGjBihXr16XfbSWtnjunfvriFDhmj9+vU6ceKEduzYoaeeekp79uypodkBvBthBrDUs88+W+4S0k033aS//vWvWr58uTp06KDf/OY3euaZZzRq1Kgqj+vv768NGzaoefPmGjhwoDp27Kg5c+Y4V3aGDBmiF198UXPnzlX79u21aNEipaamKjEx8aqez3XXXac1a9Zo165d6ty5s8aNG6exY8fqqaeeuqpxq+qZZ57RiRMn1Lp1a5d7UiRp7NixKi4u1pgxYyodIzQ0VFu3btXAgQMVHx+vp556SvPmzdOAAQMkSSNHjtT8+fO1YMECtW/fXj/84Q+VlpYm6dKKzurVq3XNNdfojjvuUN++fdWqVSutWLGi0nM6HA6tWbNGd9xxh8aMGaP4+Hjdf//9OnHihMLDw69iRgB7OExFF9QBAE5vvvmmJk6cqDNnzsjf39/T5QD4Dl9PFwAA3qqwsFDp6elKTk7Wo48+SpABvBSXmQDgMp5//nndeOONCg8P1/Tp0z1dDoDL4DITAACwGiszAADAaoQZAABgNcIMAACwGmEGAABYjTADAACsRpgBAABWI8wAAACrEWYAAIDV/h/QrjkDvPCokwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.hist(nc_scores, bins=20, range=[0, 1])\n",
    "plt.axvline(thres, color=\"black\", label=\"Threshold\")\n",
    "plt.xlabel(\"Nonconformity score\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d11fed1-c67a-4778-b993-c59224792ce6",
   "metadata": {},
   "source": [
    "#### Compute predictions sets for new data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e947b7dd-08cc-45df-82e6-b9d00c8d23d1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True,  True],\n",
       "       [ True, False, False],\n",
       "       [False, False,  True]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_nc_scores = 1 - model.predict_proba(X_test)\n",
    "\n",
    "prediction_sets = test_nc_scores <= thres\n",
    "\n",
    "prediction_sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ff023673-436a-47f7-a7d3-65fffc907c1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Versicolor' 'Virginica']\n",
      "['Setosa']\n",
      "['Virginica']\n"
     ]
    }
   ],
   "source": [
    "class_names = np.array([\"Setosa\", \"Versicolor\", \"Virginica\"])\n",
    "\n",
    "for row in prediction_sets:\n",
    "    print(class_names[row])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b492dc-92c5-48d8-8a8b-694eacf0ffee",
   "metadata": {},
   "source": [
    "# Conformal Prediction Using MAPIE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d06dd846-321a-4d04-aac4-2ad99d51a8ec",
   "metadata": {},
   "source": [
    "- MAPIE documentation: https://mapie.readthedocs.io/\n",
    "- Code is inspired by Christoph Molnar's [\"Conformal Prediction\" book](https://christophmolnar.com/books/conformal-prediction/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "763dde38-4a99-47b6-bb95-797683d1c923",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>MapieClassifier(cv=&#x27;prefit&#x27;,\n",
       "                estimator=LogisticRegression(random_state=123,\n",
       "                                             solver=&#x27;newton-cg&#x27;))</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" ><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">MapieClassifier</label><div class=\"sk-toggleable__content\"><pre>MapieClassifier(cv=&#x27;prefit&#x27;,\n",
       "                estimator=LogisticRegression(random_state=123,\n",
       "                                             solver=&#x27;newton-cg&#x27;))</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" ><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">estimator: LogisticRegression</label><div class=\"sk-toggleable__content\"><pre>LogisticRegression(random_state=123, solver=&#x27;newton-cg&#x27;)</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-3\" type=\"checkbox\" ><label for=\"sk-estimator-id-3\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LogisticRegression</label><div class=\"sk-toggleable__content\"><pre>LogisticRegression(random_state=123, solver=&#x27;newton-cg&#x27;)</pre></div></div></div></div></div></div></div></div></div></div>"
      ],
      "text/plain": [
       "MapieClassifier(cv='prefit',\n",
       "                estimator=LogisticRegression(random_state=123,\n",
       "                                             solver='newton-cg'))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mapie.classification import MapieClassifier\n",
    "\n",
    "\n",
    "cp = MapieClassifier(estimator=model, cv=\"prefit\", method=\"score\")\n",
    "cp.fit(X_calib, y_calib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6fa514ce-6bc8-477d-b241-b8c008553856",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True,  True],\n",
       "       [ True, False, False],\n",
       "       [False, False,  True]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_pred, y_set = cp.predict(X_test, alpha=1-confidence_level)\n",
    "y_set = np.squeeze(y_set)\n",
    "y_set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "acd54b34-d678-4fbb-a487-35c948eb1229",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Versicolor' 'Virginica']\n",
      "['Setosa']\n",
      "['Virginica']\n"
     ]
    }
   ],
   "source": [
    "for row in y_set:\n",
    "    print(class_names[row])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
